First- and second-order phase transitions in Ising models on small world networks, 
simulations and comparison with an effective field theory 
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We perform simulations of random Ising models denned over small-world networks and we check 
the validity and the level of approximation of a recently proposed effective field theory. Simulations 
confirm a rich scenario with the presence of multicritical points with first- or second-order phase 
transitions. In particular, for second-order phase transitions, independently of the dimension do of 
the underlying lattice, the exact predictions of the theory in the paramagnetic regions, such as the 
location of critical surfaces and correlation functions, are verified. Quite interestingly, we verify that 
the Edward- Anderson model with do = 2 is not thermodynamically stable under graph-noise. 
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I. INTRODUCTION 

Disordered systems represent one of the most impor- 
tant fields of statistical mechanics. Disorder is at the base 
of many interesting phenomena whose understanding is 
often far from being trivial or immediate. Its use varies 
from applications in condensed matter physics to com- 
puter science and, more recently, to the broad range of 
natural, artificial, and social networks. One of the most 
important analytical tool to study disordered systems is 
represented by suitable mean-field theories. Originally 
developed to understand spin glass models, the replica 
method (RM) and the cavity method (CM) are nowadays 
largely used in many other fields of statistical physics and 
computer science and represent the most powerful ana- 
lytical methods to investigate the intricate nature of the 
spin glass and other disordered phases PJ-Q. 

Among disordered models whose dimensionality - in a 
broad sense - can be considered infinite one can sin- 
gle out a "hierarchical" family of models of increasing 
difficulty such as: fully connected models (correspond- 
ing to infinite connectivity in the thermodynamic limit) 
like the Sherrington-Kirkpatrick model @ , finite connec- 

and mod- 
This lat- 



tivity models like the Viana-Bray model 
els defined over small-world networks 0, _ 
ter class of models has been introduced in recent years 
and represents an important development in modeling 
more realistic situations in which the spins, besides inter- 
acting through a random finite connectivity, distributed 
according to some given distribution with average de- 
gree c, interact also through short-range connections. In 
other words, small-world models constitute an interplay 
between purely random and regular finite-dimensional 
models. It is known that, despite the underlying finite 
dimensionality do present in these kind of graphs, in the 
thermodynamic limit, models defined on them manifest 
a mean field behavior. This fact, far from being trivial, 
to be rigorously proved, may lead to hope that the RM 
or the CM could be used to solve small-world models. 
Indeed such methods have been already successfully ex- 
ploited in [9j-fll{ for do = 1 small-world models defined 



upon adding a Poisson distributed random connectivity 
c to the underlying regular one-dimensional chain. How- 
ever, if we take a look at the mathematical structure of 
these methods we recognize the following. For what con- 
cerns the RM, we need to know analytically the two lead- 
ing eigenvectors of the transfer matrix of the Ising model 
without the shortcuts but immersed in a random external 
field; whereas for what concerns the CM, it is essential 
that the underlying graph Cq had a tree-like structure, 
i.e., no loops, at least in the thermodynamic limit. As 
a consequence, both methods seem hardly applicable to 
small- world models if do > 1- The effective field theory 
we have recently developed in [l2j], based on mapping a 
generic random model onto a non random one (the "pure 
model") SlUEl is instead applicable to these models. 
In fact, though this theory is able to give exact answers 
only in the paramagnetic (P) regions, there is no limita- 
tion in the underlying dimension do . Whereas the theory 
can be fully treated analytically for dq < 1, for do > 1, 
we can still apply it semi-analytically |3J| . All we need to 
apply the theory is to solve - analytically or numerically - 
the pure model in do dimension in the absence of the ran- 
dom shortcuts and in the presence of a uniform external 
field. The values of a certain observable Oo so obtained 
una tantum will be then used to get the corresponding 
value O for the model in the presence of the shortcuts 
and for any choice of the disorder parameters (couplings 
and connectivity). This feature, together with the fact 
that the effective field equations of the theory have a very 
simple structure and a more immediate physical interpre- 
tation compared to the equations of the RM method (in 
which the introduction of several coupled auxiliary fields 
is necessary), makes this effective field theory particu- 
larly interesting to all those applications in which do > I 
or else the number of parameters of the model is high. Of 
course one has to pay such an advantage with the impos- 
sibility to get exact results out of the P region. However, 
as we shall show in this paper, also in the other regions of 
the phase diagram the theory succeeds in giving effective 
approximations allowing us to obtain important insights 
on the frozen states, even if we do not have a direct ac- 
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cess to them, as instead the RM or the CM could do, if 
they were applicable also to models with loops. 

In this paper we consider random Ising models defined 
over small-world networks having an underlying regular 
lattice Co of dimension do = 1, 2, 3. Given the initial lat- 
tice Cq with N sites, we build the small- world network 
by adding cN/2 links uniformly spread over Co- This im- 
plies that at each site, besides the 2do neighbors, there 
are additional long-range neighbors whose number is dis- 
tributed according to a Poisson distribution with average 
c. Interactions act via a coupling Jo for the 2do short- 
range neighbors and via a further coupling J for the other 
long-range neighbors. This way of building a small- world 
network is different from the re-wiring method of Watts 
and StrogatzQ (in which the number of shortcuts per 
site is also Poissonian distributed) and it is more conve- 
nient for analytical calculations. However, just by using 
the effective field theory at the base of this paper, it is 
possible to deal with the similar re-wiring small-world 
models as well, and to show rigorously that the critical 
behavior of the two kind of models is identical [l5| . 

By using Monte Carlo (MC) simulations we check the 
predictions of the effective field theory for the critical sur- 
faces, the susceptibility, the average magnetization and 
the two point connected correlation function as a func- 
tion of the Euclidean distance r defined on Co- 

The ferromagnetic Ising model (both J and Jo posi- 
tive )on small- world networks has been extensively stud- 
ied However, of remarkable interest, for both 
its theoretical and practical implications, is the case with 
negative short-range antiferromagnetic coupling, Jo < 0. 
The anti-ferromagnetic Ising model on small world net- 
works was studied in [2l| for the case where there is only 
one antiferromagnetic coupling constant (Jo = J), but, 
apart from the fully connected case [22[ , no attention has 
been paid to the fact that, in the more general situation, 
may exist multicritical points with first- and second-order 
phase transition. In fact, when Jo < and J > 0, the 
effective field theory predicts two critical temperatures 
with first or second-order phase transitions separating 
two P regions. In this case, simulation results show large 
fluctuations and a rather slow approach to the thermo- 
dynamic limit, confirming the slow dynamics of the frus- 
trated system and the importance to have analytical or 
semi-analytical frustrated-free tools to investigate these 
models. Quite interestingly, the model with given cou- 
plings Jo and J, and slightly different values of the con- 
nectivity c, can show drastically different phase diagrams 
either showing two ferromagnetic (F) phase transitions or 
a single spin-glass (SG) phase transition. In the case of 
first-order phase transitions the magnetization disconti- 
nuity is not exactly predicted by the theory but simu- 
lation results show clearly the signature of a discontinu- 
ity for both the F and SG order parameters in corre- 
spondence of the theoretical P-SG critical temperature. 
Furthermore, we see good agreement of the susceptibil- 
ity predictions in the P phase. We also consider a two- 
dimensional modified Edwards- Anderson model [231] with 



added long-range shortcuts. This is a special case where 
besides the connectivity disorder there is also disorder 
on the value of the short-range coupling. As predicted 
by the theory, simulations confirm that any infinitesimal 
addition of shortcuts leads the system to have a finite 
temperature P-SG phase transition which coincides with 
the theoretical one. 

The paper is organized as follows. In Sees. II and III 
we recall the definition of the small- world models and the 
method, which mainly consists in finding the solution of 
the self-consistent equation (|13|) and minimizing the ef- 
fective free energy (|23|). In Sec. IV we report our MC 
simulations for several interesting cases selected for com- 
parison with the theoretical predictions. Finally, in Sec. 
V some conclusions are drawn. 



II. SMALL WORLD MODELS 

We consider random Ising models constructed by 
super-imposing random graphs with finite average con- 
nectivity onto some given lattice Co whose set of bonds 
(i,j) and dimension will be indicated by To and do, re- 
spectively. Given an Ising model - shortly the pure model 
- of N spins coupled over Co through a coupling Jo and 
with Hamiltonian 

H = - J ^2 o-i(Tj - h^o-i, (1) 

(*.i)er i 

and given an ensemble C of unconstrained random graphs 
c, c G C, whose bonds are determined by the adjacency 
matrix elements Ci.j = 0,1, we define the correspond- 
ing small-world model - shortly the random model - as 
described by the following Hamiltonian 

H c . j d = Ho - CijJijCTiO-j, (2) 

i<j 



the free energy F and the averages (O) 1 being defined in 
the usual (quenched) way as 




and (in the following a bar notation 7 indicates the two 
independent averages over the graph and couplings real- 
izations) 

W=J2 P ( C ) f^(Ui,j})(0) 1 , 1^1,2 (4) 

c6C J 

where Z c . j is the partition function of the quenched sys- 
tem 

Z c ,j = J2^ m ^ }} \ (5) 
{^} 
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(0)c;j the Boltzmann-average of the quenched system 
(note that (0)c-J depends on the given realization of the 
J's and of c: (O) = (0) c; j; for shortness we will often 
omit to write these dependencies) 



(O) 



def ^{(7;} Oc-J e 



-pH a .j{{ai}) 



(0) 



and dV ({ Ji.j}) and P(c) are two product measures given 
in terms of two normalized measures d^(Ji^) > and 
p(ci,j) > 0, respectively: 

dV{{Ji,i}) d = II Jdn(Jij) = l, (7) 



(hj),i<3 



p{c) = n y, pfe)=i- (8) 

{i,j),i<J c il3 -=0,l 

The variables Cjj G {0, 1} specify whether a "long-range" 
bond between the sites i and J is present (cjj = 1) or 
absent (cjj = 0), whereas the Jij 's are the random vari- 
ables of the given bond For the Cjj's, we shall 
consider the following distribution 



P( c ij) = ^c„,i + (l - ^) 5 c „,o, 



(9) 



where c > 0. This choice leads in the thermodynamic 
limit N — > oo to a number of long range connections per 
site distributed according to a Poisson law with mean 
connectivity c. 

In this paper for the Jjj 's we will assume the distribu- 
tion 



dJi j 



S (Jij - J) , 



(10) 



For the short-range nearest-neighbor coupling, J we 
will consider the distribution, 



dJ 



5(J - a), 



(11) 



except in the last case studied, the modified Edwards- 
Anderson model, where we consider, 

^ 1 = ^o-a) + ^(J + a) . (12) 



III. AN EFFECTIVE FIELD THEORY 

Depending on the temperature T, and on the param- 
eters of the probability distributions, d[i and p{ci,j), the 
random model may stably stay either in the paramag- 
netic (P), in the ferromagnetic (F), or in the spin glass 
(SG) phase. In our approach for the F and SG phases 
there are two natural order parameters that will be indi- 
cated by and m^- SG \ Similarly, for any correlation 



function, quadratic or not, there are two natural quanti- 
ties indicated by and C^- SG \ and that in turn will 
be calculated in terms of and m^ SG \ respectively. 
To avoid confusion, it should be kept in mind that in our 
approach, for any observable O there are - in principle - 
always two solutions that we label as F and SG, but, for 
any temperature, only one of the two solutions is stable 
and useful in the thermodynamic limit. 

In the following, we will use the label o to specify that 
we are referring to the pure model with Hamiltonian (TTJ . 
Let mo((3Jo, f3h) be the stable magnetization of the pure 
model with coupling Jo and in the presence of a uni- 
form external field h at inverse temperature /?. Then, 
the order parameters m^, E=F,SG, satisfy the follow- 
ing self-consistent decoupled equations 



m (E) = m ( / 3J ( f ) ,/3,/ (s) m (s) + /3/i), 



(13) 



where the effective couplings j( F ) , j( SG ) , and Jq SG ^ 
are given by 



p.J^ = c / d/*(J w ) tanhGSJy), 



/3J( SG ) = c / ^(J^tanh^Jij). 



(14) 



(15) 



For a constant short-range coupling distributed as in fp~T|) 



4 F) = « 

Pjf G) = tanh^tanh 2 ^)). 
and for the bimodal distribution 03 



(16) 



J, 



(F) 







/?J SG = tanh _1 (tanh 2 (/3a)) 



(17) 



For the correlation functions we have C^, £=F,SG, 
where 

= C {fiJ?\pj<n m ™ +ph) + (i) , (18) 

where Cq((3Jq, fih) is the correlation function of the pure 
model. For the corrective Q(l /N) term in Eq. (fT5f we 
remind the reader to Eq. (33) of [12] . Let us indicate by 
C« and the averages and the quadratic averages 
over the disorder of the correlation function of degree, 
say k. Then, C (1) and C (2) , are related to C (F) and 
C (SG \ as follows 

C (i) = cf (F) j inF; (19) 

C (1) = 0, kodd, inSG, (20) 
C (1) = C (SG) , /ceven, in SG, (21) 
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and 



^(SG) 



in F, 



in SG. 



(22) 
(23) 



In particular, for the susceptibility x^ of the random 
model we have: 



1 



Xo (PJP, 0J (E) m< E > + ph 



(24) 



where xo stands for the susceptibility xo 01 the pure 
model divided by p (we will adopt throughout this di- 
mensionless definition of the susceptibility) and similarly 
for the random model. For the case £ =F without disor- 
der (dfi(J') = 8{J' - J)dJ' and d^ (J ) = S(J - a)dJ ), 
Eq. ([24]) was already derived in |2(| by series expansion 
techniques at zero field (h — 0) in the P region (where 
to = 0). 

Among all the possible stable solutions of Eqs. (fT3")) . in 
the thermodynamic limit, for both E=F and £=SG, the 
true solution m' s ), or leading solution, is the one that 
minimizes L^' where 



def /?J (S) {mf 



■ Ph{pJ^\pJ^ ) m + ph)^b) 



fo (PJo j fiti) being the free energy density in the thermo- 
dynamic limit of the pure model with coupling Jo and in 
the presence of an external field h, at inverse tempera- 
ture p. A necessary condition for a solution to be 
the leading solution is the stability condition: 

X0 (pWjP,pjWmW+ph)pWjW < 1. (26) 

For the localization and the reciprocal stability be- 
tween the F and SG phases we remind the reader to Sec. 
HID of [IH . We recall however that, at least for lattices 
£ having only loops of even length, the stable P region 
is always that corresponding to a P-F phase diagram, so 
that in the P region the correlation functions must be 
calculated only through Eqs. (fH?|) and (|2"2"j) . 



exact also in the limit c — > + , in the case of second- 
order phase transitions, and in the limit c — > oo (see Sec. 
IIIC of [13] )• Note also that the order parameters m^ s \ 
and then the correlation functions, are by construction 
always exact in the zero temperature limit. 



IV. SIMULATIONS AND COMPARISON WITH 
THE THEORY FOR GIVEN COUPLINGS 

The Monte-Carlo simulations presented in this work 
were made using a local spin-flip dynamics with a 
Metropolis acceptance probability |24j |. 

Throughout this work we estimate the susceptibility, 
in the P phase by x. = N(m 2 ), and in the ferromagnetic 

phase by x = N ((m 2 ) - (|m|) 2 j , where m = YaLi a i 

is the magnetization of the system and N — L d ° is the 
total number of spins in the lattice of side L. The Binder 
cumulant[25], defined by 



U L = 1- 



(to 4 ) 



(28) 



was used to locate the critical points. The cumulants 
Ul and Ul 1 , for two systems of different sides L and L', 
plotted as a function of temperature, cross at the critical 
point at a value, U* that characterizes the universality 
class of the model. 

To study spin-glass phases we calculate the overlap 
order-parameter, q — YiLi °i '^i^ obtained from two 
replicas of the system with spins and of' 1 '. The ob- 
served distribution of the values of q is measured for a 
given realization of the disorder which corresponds to 
taking a thermal average. Subsequently, by considering 
different samples, an average over disorder is done: 



P{q) 



i M 

3 = 1 



(29) 



where qj is the value of the overlap parameter at time 



The inverse critical temperature P, 
following exact equation 



Xo (pMjf\0)pWjW = l, pW<P%\ (27) 



is solution of the s ^ e P b ^ s the Kronecker delta, M is the number of 
simulation Monte Carlo steps (MCS) after thermal equi- 
libration is reached, and the bar denotes averaging over 
disorder. The Binder cumulant for the overlap order pa- 
rameter can be defined by 



where fi cQ is the inverse critical temperature of the pure 

model with coupling Jq . When Jo > 0, the constrain in 
Eq. (|27p ensures the uniqueness of the solution. However, 
if Jo < 0, Eq. (|2T1) in general admits either or at least 
2 solutions (in principle also 4, 6, etc.). 

We end this section by stressing that this method is 
exact in all the P region and, at least for second-order 
phase transitions, provides the exact critical surface, be- 
havior and percolation threshold, and that, in the ab- 
sence of frustration, the order parameters become 



U q ,L = 1 



<9 4 ) 



(V) 



(30) 



We study the small-world model with the distribution 
of random bonds defined in Eq. © and a fixed positive 
long-range coupling constant as in Eq. (|T0|) . We start 
in subsection IIV Al to study the ferromagnetic case with 
Jo > and the location of the critical points for differ- 
ent values of c and do = 1,2 and 3. In subsection IIV Bl 
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we analyze the spin-spin correlation function above the 
critical temperature for d = 2. In subsection IIV CI for 
the special case Jo = we study the magnetization and 
susceptibility at non-zero external field above the criti- 
cal temperature. In subsection IIVDI we consider a one 
dimensional system with Jo negative where two second- 
order phase transitions are predicted. In subsection IIV El 
we analyze the same one-dimensional model for couplings 
and connectivity such that either two first-order phase 
transitions or a spin-glass phase are predicted. Finally 
in subsection IIV Fl we study a two dimensional Edwards- 
Anderson model with added long-range shortcuts. 



A. The Paramagnetic-Ferromagnetic line of critical 

points 

From the susceptibility of the pure Ising model in a 
hypercubic lattice of dimensionality do and Eq. (|2"7|) we 
obtain the location of the P-F line of critical points in 
the c-T plane. For do = 1 we use the known analytical 
expression for the susceptibility, xo ((3Jo, 0) = exp(2/3 J ), 
that applied to Eq. (f2"Tl) reproduces the same formula of 
[j| for the P-F and P-SG lines. For higher dimensions we 
use numerical results obtained from Metropolis Monte- 
Carlo simulations. 

For do = 2 the pure model susceptibility was de- 
termined for several systems sides up to L — 128 to 
check for finite-size effects. The pure model was sim- 
ulated for 40 temperatures, in the range 0.1 < /3Jo <= 
/3 c o Jo = 0.44068.... For d = 3 we studied systems of side 
L = 8 and L = 16 also for 40 temperatures in the range 
0.05 < f3J <= /3 c0 Jo = 0.2216546(10) (H- In all of the 
simulations reported we neglected the first 10 5 MCS/N 
and made measurements in the remaining 10 6 MCS/N 
steps. In Fig. Q]we plot the critical lines for do = 1 with 
J I Jo = 3/5, and for do = 2 and 3 with J/ Jo = 1. Note 
that in the figure, for do = 2 and 3, there are several lines 
corresponding to the use of susceptibility estimates ob- 
tained from systems of different size. Nevertheless, in the 
scale of the plot no finite size effects can be seen. This is 
a consequence of the fact that only very near c = the 
solution of Eq. (|2T|) uses values of the pure model suscep- 
tibility near the critical point. Far from the critical point 
of the pure model the susceptibility and consequently the 
estimates of the critical line in the disordered model do 
not show finite-size effects. 

In order to compare the predictions based on Eq. (|27[) 
we measured by direct simulation the location of the crit- 
ical points for several values of the average connectivity, 
c = 0, 0.5, 1, 1.5, 2, 5 and 10 and the results, for each spa- 
tial dimension, are plotted in Fig. [TJ The average over 
disorder was done by considering averages over 10 sam- 
ples. 

For the case do = 2 we present in Table fl] detailed nu- 
merical results. For the cumulant crossing Pl.L'Jq (third 
column) listed for several values of c we estimate a statis- 
tical error 0.0005 which allows us to claim a good agree- 




FIG. 1: (Color online) Lines of critical points in the /3 — c 
plane for do = 1 (O) with J / Jo = 3/5, and for do = 2 (□) and 
do = 3 (*) with J I Jo = 1 as obtained from the pure model 
susceptibility and Eq. (|27|) . The symbols are estimates of 
the critical points obtained from simulation and by using the 
Binder cumulant intersection technique. 



mcnt between the simulations and the theoretical predic- 
tion obtained from Eq. (|2"T|) (fourth column). The values 
of the cumulant at the critical point for the pure model 
(c = 0) are close to the value 0.61069 calculated in [27| . 
The long-range links introduced by the disorder change 
the universality class from 2d Ising to mean-field. 

The mean-field value of the Binder cumulant at criti- 
cality for an infinite system is predicted to be 0.2705[28l 
[29|,[3l[ which is close to our estimates of Binder cumulant 
intersections listed in table fl] (fifth column) for which we 
estimate an error equal to 0.02. In three dimensions we 
simulated only systems of side L = 8 and 16. For c = 
the two cumulants intersect at the value Ul = 0.486 close 
to the estimation 0.46521 reported in [26]. For the other 
values of c studied, c = 0.5, 1, 1.5,2,5 and 10 the inter- 
section was measured near Ul = 0.31. In one dimension 
we studied only c = 0.5, 5 and c = 10 and also the Binder 
cumulant intersections were found to be near 0.3 for in- 
tersections of L = 128, 256, 512, 1024 with V = 2048. 

Furthermore, for do = 2, we measured the scaling with 
system size of the a verage value of the absolute value 
of the magnetization (\m\)(/3 c Jo) and the susceptibility, 
xiftcJo) at the critical point. For a mean-field universal- 
ity class these quantities are expected to scale at critical- 
ity like (jmj)^ ~ i\r 1/4 _ L -d /i and ^ c _ N i/2 ^ L d /2 
where N is the total number of spins 0, [29|, [3l|. In 
Fig. [2] we show these quantities plotted in bi- logarithmic 
scale as a function of system size for the different c values 
studied. For the magnetization, the slope of the straight 
line fit for c = is —0.125(1) and for the susceptibility 
is 1.752(1), consistent with the known exact exponents 
of the pure do=2 Ising model. For c = 0.5,1,1.5 and 
2 we got, respectively, for the magnetization exponents 
-0.49(3), -0.48(3), -0.52(2) and -0.48(3) close to the 
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TABLE I: (Color online) Results of Binder cumulant cross- 
ings for do = 2 and Jo > 0. In the column Pl,l'Jq we list 
the crossing inverse temperatures between the cumulants for 
systems of side L with L' = 128. In the column j3 c Jo we list 
the critical parameters obtained from Eq. (I27|l . The crossing 
values of the cumulant are listed in the column (72 L , . 



c 


L 


Pl.l'Jo 




Ul, 128 




16 


0.4410 




0.612 





32 


0.4409 


0.440686... 


0.612 




64 


0.4411 




0.615 




16 


0.2968 




0.33 


0.5 


32 


0.2964 


0.2963 


0.30 




64 


0.2960 




0.27 




16 


0.2466 




0.32 


1.0 


32 


0.2461 


0.2461 


0.28 




64 


0.2460 




0.28 




16 


0.2139 




0.30 


1.5 


32 


0.2137 


0.2134 


0.29 




64 


0.2136 




0.28 




16 


0.1894 




0.30 


2.0 


32 


0.1891 


0.1897 


0.27 




64 


0.1888 




0.24 




16 


0.1173 




0.28 


5.0 


32 


0.1174 


0.1167 


0.29 




64 


0.1173 




0.27 




16 


0.0730 




0.28 


10.0 


32 


0.0730 


0.0728 


0.28 




64 


0.0731 




0.30 


expected value 


—0.5; whereas for the susceptibility we ob- 



tain the exponents 1.03(5), 1.06(7), 0.97(2) and 1.01(5), 
also close to the expected result for mean-field behavior. 



B. Correlation functions above the critical 
temperature 

From Eq. (fTS]) we see that the effective field theory, 
in the P region and zero external field, predicts the spin- 
spin correlation function of the random model to 
be, in the thermodynamic limit, equal to the correlation 
function of the pure model calculated at the same tem- 
perature. In order to check this result we calculated, 
from simulation, 

C( F )(r) = (^y, (31) 

where uq is an arbitrary spin and ay is one spin at 
Euclidean distance r from the spin <7o, measured on 
the lattice Co. We considered the case do = 2 with 
J = Jo and c = 1 and 2 and the inverse temperatures 
/3Jo = 0.1,0.17571. These temperatures are above the 
critical temperatures for the two values of c studied. We 
studied also the correlation function at the critical in- 
verse temperatures /3 c Jo = 0.2461 and /3 C Jo = 0.1897 for 
c = 1 and c = 2, respectively (see Table H]). These calcu- 
lations were done for system sides L = 8, 16, 32 and 64. 
In Fig. |3] we can see that, as the system size increases, 
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FIG. 2: (Color online) In (a) we plot the susceptibility, \ c , 
and in (b) the magnetization, {\m\} , at the critical point 
for the do = 2 model with J/ Jo = 1 as a function of the 
system side L. Both in (a) and (b) the data are for c = 
0(*),0.5(x),l(O),1.5(+),2(D). 

the curves C F (r), for c > 0, approach the data points for 
c = calculated for a system of size L = 64. Note that, 
as we approach the critical temperature, the correlation 
function at large r reaches a finite value that decreases 
as the system size increases. This finite constant is just 
due to the finite size of the system and it is predicted by 
the theory to vanish as 1/N, in the thermodynamic limit, 
but at the same time it is responsible for the divergence, 
with system size, of the susceptibility at the critical point 
(see Eq. (33) and subsequent comments of Ref. fl2j]). 

C. Susceptibility and Magnetization at non-zero 
external field above the critical temperature 

In the P phase at zero external magnetic field the effec- 
tive field theory prediction for the susceptibility is exact, 
consistently with the exact predictions for the critical 
temperatures. The question remains whether the sus- 
ceptibility prediction is a good approximation for non- 
zero field. To verify this we made simulations for the 
case Jo = with a positive long-range coupling (as in 
Eq. (|10p) (in other words the simplest version of Viana- 
Bray model). For this particular case the magnetization 
is predicted to be given by the solution of the following 
equation 

to = tanh (c to tanh(/3J) + @h) (32) 
and the susceptibility (divided by /3) is given by, 

1 — TV? 

X(^^)= 1 _ c(1 _ m2)tanh ^ (33) 
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FIG. 3: (Color online) Correlation functions C' F ^(r) as a 
function of Euclidean distance r for J = Jo and c = 1 (in 
(a), (b) and (c)) and c = 2 (in (d), (e) and (f)) at different 
temperatures, /3J = 0.1 ((a) and (d)), 0.17571 ( (b) and (e)). 
In (c) with /3 C Jo = 0.2461 and in (f) with /3 C J = 0.1897 (the 
critical temperatures for c = 1 and c = 2, respectively). The 
lower line (which is blue in the online version of this paper) 
is C {F \r) for the pure model c = and for a system side 
L = 64 (0). Data for four system sides are plotted, L = 8 
(O), L = 16 (+), L = 32 (x), L = 64 (0). 



In Fig. |4]we compare the results of the above predic- 
tions for the magnetization and susceptibility with sim- 
ulation results for c = 2. For averaging purposes we con- 
sidered 10 samples. The plots correspond to three tem- 
peratures (/3J)- 1 = 1.2903 (first row), ((3J)- 1 = 1.8182 
(middle) and {(3J)' 1 = 2.1739 (bottom). The critical 
temperature is (/3 C J) _1 = 1.8205. By construction, in 
the limit of strong field the magnetization prediction be- 
comes exact and similarly in the limit of small field above 
the critical temperature. In the intermediate field range 
we see that the magnetization prediction and simulation 
results in general do not agree. However, above the criti- 
cal temperature, the susceptibility obtained from simula- 
tion and the theoretical prediction given by Eq. ([33)) are 
very close to each other over the full range of field val- 
ues studied. Note that for (/3J) _1 = 1.8182, close to the 
critical temperature, the simulation susceptibility shows, 
as expected, a strong finite-size effect at zero field. 



FIG. 4: (Color online) Magnetization ( (a), (b) and (c) ) 
and Susceptibility ( (d), (e), (f)) for the Viana-Bray model 
and c = 2 for three temperatures as a function of the ex- 
ternal magnetic field. The lines are the theoretical predic- 
tions ( see Eqs. (|32p and J33}) and the data points are sim- 
ulation results for N = 128(.), 256(o), 512(+), 1024(x) and 
2048(0)- For the case studied here the critical temperature 
is, (/3 C J) _1 = 1.8205. 



It is worth to observe that with respect to our effective 
field theory, the Viana-Bray model represents the worst, 
i.e., the most difficult, case. The theory in fact, by con- 
struction, takes exactly into account all the effects due to 
the short-range couplings and to the short-loops present 
in the given lattice Cq, and the greater is do, the greater 
is the level of accuracy of the theory also out of the P 
region (at least in the absence of frustration), while in 
the Viana-Bray model topologically we have do = 0. 



D. Negative short-range coupling and second order 
phase transitions 

In the case Jo < the theory allows for the occur- 
rence of two second-order phase transitions. At low 
and high temperatures the system is disordered and 
in the intermediate temperatures a ferromagnetic phase 
arises. The simulations confirm this phase diagram pic- 
ture. We made simulations at do = 1) for Jo = —0.5, 
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FIG. 5: (Color online) Magnetization (top) and susceptibility 
(bottom) as a function of temperature for a do = 1 system 
with Jo < 0. The average connectivity is c = 1.4 and Jo = 
—0.5, J/ Jo = —20. The line is the theoretical prediction 
and the data points correspond to different system sizes, L — 
512(.), 1024(o), 2048(+), 4096(x), 8192(o), 16384(A). 



J/ Jo — —20 and c = 1.4. We made averages over dis- 
order by considering 50 samples and we studied system 
sizes L = 512,1024,2048,4096,8192 and 16384. In Fig, 
[5] we plot the magnetization and the susceptibility as a 
function of temperature, T. The two critical points are 
predicted to occur at T cA = 2.985 and T c , 2 = 9.207. The 
predicted values of the magnetization in the intermediate 
temperature range, T Cj i < T < T c ^ 1 are different from 
the simulation results (see Fig. [3]). However, the theo- 
retical predicted susceptibility, in both P phases, is very 
close to the simulation susceptibility approaching each 
other as the system size increases. 

Note that this model is a frustrated system so that 
large fluctuations and strong finite size effects are present, 
especially close to the lower temperature critical point 
where we do not perform high precision simulation as 
it requires averaging over a large number of samples. 
We have studied, in detail, the high temperature criti- 
cal point where we applied the cumulant crossing tech- 
nique. The intersection temperatures of the cumulants 
for L = 512,1024,2048,4096 with the system of size 
L = 16384 were, 9.24,9.35,9.23,9.32, respectively, from 
which we can estimate a critical temperature equal to 
9.29(6). This estimate is close but slightly higher than 
the predicted critical temperature, T c ^ = 9.207. How- 
ever, considering the statistical error and the finite size 
corrections we cannot exclude a convergence toward the 
theoretical value in the thermodynamic limit. The cor- 
responding intersection values of the cumulant were, 
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FIG. 6: (Color online) Binder cumulant Ul as a function of 
temperature for a do = 1 system with Jo < 0. The average 
connectivity is c = 1.4 and Jo = —0.5, J/ Jo = —20. The ver- 
tical line is the theoretical prediction for the critical tempera- 
ture and the data points correspond to different system sizes, 
L = 512(0, 1024(c), 2048(+), 4096(x), 8192(o), 16384(A). 



0.23,0.21,0.23,0.21. These values are smaller than the 
values in the range 0.27 — 0.3 that we measured in sec- 
tion EES for J > 0. In the Fig. ©we plot the Binder 
cumulant as a function of temperature for the system 
sizes studied. 



E. Negative short-range coupling, first-order and 
spin-glass phase transitions 

We considered the case for c = 10, Jo = —0.9 and 
J = 0.5 at do = 1. From the theory it turns out that 
for temperatures above 2.38 only the zero magnetization 
solution is stable but for temperatures lower than this 
value there are always two stable solutions, one with non- 
zero magnetization and another with zero magnetization. 
For temperatures T < 2.34 the nonzero magnetization 
solution has the lower free-energy so that a first-order 
phase transition is predicted at this temperature. The 
theory also predicts a possible spin-glass phase transition 
at a temperature, T Cj sg = 1-88. Note that the theory 
always predicts continuous spin-glass phase transitions. 

We made simulations for systems of size L — 512, 
1024, 2048, 4096 and 8192 and as before we neglected 
10 5 MCS/N for equilibration purposes and we made mea- 
surements for 10 6 MCS/N. The averaging over disorder 
was made by considering 50 samples. The results show 
a first order phase transition occurring at slightly lower 
temperatures than the one predicted by the theory. The 
probability distribution of the magnetization clearly ex- 
hibits (see Fig. [7]) the behavior characteristic of first- 
order phase transitions (3fj| namely the emergence, at 
temperatures close to the transition, of two maximum lo- 
cated at symmetric nonzero values together with a third 
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FIG. 7: (Color online) Simulation results for the magneti- 
zation probability distribution, P(m), for the model with 
c = 10, Jo = —0.9 and J = 0.5 at do = 1. In (a) we plot 
curves for L = 8192 (solid line), L = 4096 (dotted line) at 
T = 2.05 and L = 2048 (dashed line) at T = 2 and in (b) we 
plot the corresponding data for L = 8192 and L = 4096 at 
T = 1.86 and L = 2048 at T = 1.83. 



maximum near zero magnetization. 

In Fig. [UJa) we show the simulation results for the av- 
erage magnetization, the magnetic susceptibility (b), and 
the magnetization Binder cumulant (c). The simulation 
susceptibility follows very closely the zero magnetization 
theoretical susceptibility up to temperatures lower than 
the predicted T c . Note that, in the case under study, 
the zero magnetization solution is a stable solution at 
any temperature and the theoretical prediction of the 
phase transition is based on comparison of the value of 
free-energies of the solutions. The theory does not pre- 
dict correctly the location of the transition since values 
of the free-energy of the non-zero magnetization solu- 
tion are not given exactly by the theory. The simulation 
Binder cumulant shows, near the transition temperature, 
the expected increasingly negative values, as the system 
size increases [3(| • Interestingly, the simulation data give 
a transition temperature close to the theoretically pre- 
dicted spin-glass transition temperature. 

We also studied the overlap order parameter distribu- 
tion and the results are shown in Fig. [SJ Here, also, the 



FIG. 8: (Color online) Simulation results for the average 
magnetization (a), magnetic susceptibility (b) and Binder cu- 
mulant Ul (c) as a function of temperature for the model 
with c = 10, Jo = —0.9 and J = 0.5 at do = 1. In 
(b) we have used x — N< m 2 > at any temperature. The 
vertical lines are at T c = 2.34 and T c ,sg = 1.88. In (b) 
the line is the theoretical susceptibility for the zero magne- 
tization solution. The simulation data points are for L = 
512(0), 1024(x), 2048(.), 4096(+), 8192(o). 



behavior expected for a first-order transition temperature 
was observed. The negative value minima of the Binder 
cumulant for the overlap order parameter are steeper and 
occur at slightly higher temperatures, for the same sys- 
tem sizes, as compared with the corresponding quantity 
for the magnetization. 

Quite interestingly for c = 4.5 and the values of the 
coupling constants, J = —0.9 and J = 1, the theory 
predicts the absence of ferromagnetic phase transitions 
and only a continuous P-SG phase transition located at 
T C (5G) = 2.29. 

The obtained simulation results are shown in Fig. 1101 
The number of samples was 100 and the simulation times 
considered here were the same as for the case c = 10. In 
the top plot we see that U q x for different system sizes 
intersects very near the theoretically predicted critical 
temperature. The inset of the top plot shows the average 
overlap order parameter that decreases with system size 
above the spin-glass critical temperature. In the lower 
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FIG. 9: (Color online) Simulation results for the average 
overlap order parameter, < \q\ > (a), \sg = N< q 2 > (b) 
and Binder cumulant, U q (c), as a function of tempera- 
ture for the model with c = 10, Jo = —0.9 and J = 0.5 
at do = 1. The data points are simulations for L = 
512(0), 1024(x), 2048(.), 4096(+), 8192(o). 
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FIG. 10: (Color online) Results for the model with c = 4.5, 
Jo = —0.9 and J — 1 at do = 1- In the top plot we show 
the temperature dependence of the Binder cumulant for the 
overlap order-parameter. The lines are fits made for each sys- 
tem size. The vertical line is the location of the predicted 
spin-glass critical temperature. The inset shows the aver- 
age values of the overlap order parameter together with the 
effective field theory prediction for this quantity. In the bot- 
tom plot we show the temperature dependence of the mag- 
netic susceptibility. The inset shows the average magneti- 
zation as a function of temperature. For all the plots the 
data points correspond to simulations for systems of sides, 
L = 512(0), 1024(x), 2048Q, 4096(+). 



plot we see that the magnetic susceptibility is almost 
independent of system size and it agrees with the theo- 
retical prediction for temperatures above the spin-glass 
phase transition while it deviates at lower temperatures. 

F. Bimodal Edwards-Anderson model 

The bimodal Edwards- Anderson model 23] is a disor- 
dered spin model where the nearest neighbor coupling Jo 
of Ising type spins has the bimodal distribution in (1121) . 
Here, we consider the two dimensional square lattice ver- 
sion of the model with additional long-range couplings. 
The pure model (without long range shortcuts) is known 
to show a P-SG phase transition only at zero tempera- 
ture being the lower critical dimension of the model equal 
to two [32], [33[ . 

To apply the effective field theory, for each phase E =F 
or SG, we have to consider the pure Ising model magneti- 
zation (|13p and susceptibility (|2"4")l calculated by using the 
definitions of the effective long- and short-range couplings 
j( s ) and jff\ from Eqs. ([T4 ) -([P7 |) . respectively. By us- 



ing the numerical two-dimensional Ising model magnetic 
susceptibility, we can obtain the expected location of the 
spin-glass phase transition. The result of this calculation 
is shown in Fig. QTJ We stress that theory predicts that 
the inclusion of an arbitrary small number of long-range 
shortcuts in the Edwards- Anderson model leads always 
to a finite temperature phase transition. In fact, in the 
limit of an infinitesimal addition of short-cuts, the the- 
ory predicts a P-SG transition at the finite value given 
by (3 c a = tanh -1 [ v /tanh(O.44O68..0] = 0.7642..., where 
0.44068... is the critical inverse temperature of the regu- 
lar two-dimensional Ising model with a unitary positive 
coupling. This implies that the EA model with short- 
cuts, in the limit c — > + , is not equivalent to the orig- 
inal EA model without short-cuts. In other words, the 
EA model is not thermodynamically stable under graph 
noise. 

We made simulations for the EA model with a = 1, 
J = 1 and c = 1 to compare with the predictions 
of the theory. The average over disorder was done by 
considering 100 samples. We neglected the first 10 5 
MCS/N for equilibration and made measurements for 10 6 
MCS/N. For c = 1 the expected critical temperature is 
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FIG. 11: (Color online) Phase diagram for the do — 2 
Edwards- Anderson model with long-range connections and 
coupling constant J = 1 . The curve is the value of f3 c a for 
the SG phase transition as obtained from equation (|26[) and 
our numerical results for the pure Ising model susceptibility 
for a system of side L — 128. 

T ( C SG) = 1.6057. In Fig. [12] (top panel) we obtained 
crossings of the overlap parameter Binder cumulant for 
system sizes L — 8, 16 and 32. The statistical error of the 
Binder cumulant, as expected, increases with system size 
and we excluded from the Binder cumulant intersection 
calculations the data for L = 64. Our numerical esti- 
mate of the critical temperature is T c = 1-6(1) still 
consistent with the expected critical temperature. 

Since, for an Ising model with zero coupling and 
zero external magnetic field xo = 1, the prediction 
for the magnetization and the magnetic susceptibility is 
< | m | >= and x = [1 - ctanh^J)]" 1 . In Fig. QJ] 
(lower panel) we plot the magnetic susceptibility simu- 
lation results together with this theoretical prediction. 
We see that in the P phase the numerical estimates of 
the susceptibility approach the theoretical curve as the 
system size increases. 



V. CONCLUSIONS 

In this work we have compared the predictions of an 
effective field theory for several Ising models on small 
world networks with Monte Carlo simulation results. All 
the predictions of the theory, where it is known to be 
exact, namely in the P region at zero external field, were 
confirmed by the simulation results. In particular we 
have checked the critical surfaces (T, c), the temperature 
dependence of the susceptibility, and the spin-spin corre- 
lation function. 

Furthermore, for the simplest version of the model, 
i.e., the Viana-Bray model (where do = 0), we stud- 
ied the effect of a non-zero external field. Although the 
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FIG. 12: (Color online) Results for the do = 2 Edwards- 
Anderson model with long-range connections with a = J = 1 
at c = 1. The top plot shows values of the Binder cumulant 
as a function of temperature and the bottom plot shows the 
temperature dependence of the magnetic susceptibility. The 
simulations were done for L = 8 (O)i L — 16 (x), L = 32 
(.) and L = 64 (+). In the top plot the vertical line is the 
theoretical prediction of the spin-glass critical temperature 
and the lines are polynomial fits of the Binder cumulants for 
each system size. In the lower plot the line is the theoretical 
prediction for the magnetic susceptibility. 



theory is not exact in this case, there is a reasonable 
agreement between the predicted field dependence of the 
magnetization and susceptibility specially above the crit- 
ical temperature. For the case of J/ Jo = —20 and c = 1.4 
at do = 1 we verified the existence of two second order 
phase transitions. A good agreement between the the- 
oretical temperature dependence of the susceptibility in 
the P phases and simulation results was observed. The 
location of the high temperature critical point was ex- 
plicitly verified using the cumulant crossing technique. 

For a one-dimensional model with Jo = —0.9, J = 0.5 
and c = 10, we observed a first order phase transition 
as predicted by the theory but at a lower temperature. 
This may be a consequence of the fact that in this case, 
unlike the second-order phase transitions cases, the zero 
magnetization solution does not become unstable at the 
transition temperature and - consistently with the fact 
that the theory does not give exact results out of the 
pure P regions - the critical point obtained as the point 
where the two free-energy values equal is not exactly pre- 
dicted. Quite interestingly, we find that also the P-SG 
phase transition turns out to be first-order and, further- 
more, its critical point seems to coincide with the theo- 
retical one. 
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With the couplings Jo = —0.9, J = 1 and a smaller av- 
erage connectivity, c = 4.5 the theory predicts that only 
a P-SG transition is present. We studied by simulation 
this spin-glass critical behavior and we found the critical 
temperature very close to the theoretical prediction. 

We also introduced a model not previously studied, the 
two dimensional Edwards- Anderson model with added 
long-range shortcuts, where we confirmed that even an 
infinitesimal inclusion of shortcuts makes the spin-glass 
phase transition to occur at a hnitc non-zero tempera- 
ture. In other words, as the theory predicts, we find 
that the two dimensional Edwards- Anderson model is not 
thermodynamically stable under graph-noise. 

The class of disordered models for which the theory 
is applicable is very wide and its application just relies 
on the availability of numerical or analytical results for 
the susceptibility of non disordered models in arbitrary 



do dimensions. When there is strong frustration, sim- 
ulations are difficult to perform requiring large simula- 
tion times and averages over many samples. Our results 
clearly confirm the usefulness of the effective field the- 
ory proposed in [l2j by giving accurate predictions for 
the models phase diagram. The possibility to improve 
the theory out of the P region opens a new interesting 
challenge. 
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